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In  this  paper,  we  construct  high-order  hyperbolic  residual-distribution  schemes  for  gen¬ 
eral  advection-diffusion  problems  on  arbitrary  triangular  grids.  We  demonstrate  that  the 
second-order  accuracy  of  the  hyperbolic  schemes  can  be  greatly  improved  by  requiring  the 
scheme  to  preserve  exact  quadratic  solutions.  We  also  show  that  the  improved  second- 
order  scheme  can  be  easily  extended  to  the  third-order  by  further  requiring  the  exact¬ 
ness  for  cubic  solutions.  We  construct  these  schemes  based  on  the  Low-Diffusion-A  and 
the  Streamwise-Upwind-Petrov-Galerkin  methodology  formulated  in  the  framework  of  the 
residual-distribution  method.  For  both  second-  and  third-order-schemes,  we  construct 
a  fully  implicit  solver  by  the  exact  residual  Jacobian  of  the  second-order  scheme,  and 
demonstrate  rapid  convergence  of  10  15  iterations  to  reduce  the  residuals  by  10  orders 
of  magnitude.  We  also  demonstrate  that  these  schemes  can  be  constructed  based  on  a 
separate  treatment  of  the  advective  and  diffusive  terms,  which  paves  the  way  for  the  con¬ 
struction  of  hyperbolic  residual-distribution  schemes  for  the  compressible  Navier-Stokes 
equations.  Numerical  results  show  that  these  schemes  produce  exceptionally  accurate  and 
smooth  solution  gradients  on  highly  skewed  and  anisotropic  triangular  grids,  including 
curved  boundary  problems,  using  linear  elements.  We  also  present  Fourier  analysis  per¬ 
formed  on  the  constructed  linear  system  and  show  that  an  under-relaxation  parameter  is 
needed  for  stabilization  of  Gauss-Seidel  relaxation. 


I.  Introduction 

In  many  flow  simulations,  accurate  predictions  of  solution  gradients,  such  as  velocity  and  temperature 
gradients,  are  essential  for  design  and  analysis  purposes  as  they  are  directly  related  to  the  physical  quantities 
of  interest:  the  viscous  stresses,  the  vorticity,  and  the  heat  fluxes.  However,  it  is  widely  accepted  that  accurate 
and  smooth  solution  gradients  cannot  be  achieved  with  conventional  schemes  on  fully  irregular  unstructured 
grids.1,2  In  conventional  schemes,  the  gradients  are  obtained  typically  with  a  lower  order  of  accuracy 
(typically  through  reconstruction),  and  they  are  also  subject  to  numerical  oscillations  on  such  grids.  The 
resolution  of  this  issue  is  very  important  for  justifying  the  use  of  high-fidelity  models  in  engineering  design, 
analysis,  and  optimization,  especially  for  applications  involving  complex  geometries.  The  ability  to  predict 
the  gradients  on  irregular  grids  is  even  more  critical  for  grid  adaptation,  a  vital  technique  for  efficient  CFD 
calculations  in  high-order  methods,3  because  the  grid  adaptation  almost  necessarily  introduces  irregularity 
in  the  grid.  In  fact,  current  practices  in  grid  adaptation  often  avoid  adaptation  in  certain  regions  such  as 
within  boundary  layers  where  grid  irregularity  has  a  severe  impact  on  the  solution  quality.4  Therefore, 
numerical  schemes  that  can  accurately  predict  solution  gradients  on  irregular  grids  need  to  be  developed,  so 
that  the  power  of  grid  adaptation  can  be  fully  exploited. 

The  first-order  hyperbolic  system  method,  or  the  hyperbolic  method  for  short,  which  was  proposed  in  2007 
for  diffusion,5  enables  the  construction  of  such  schemes.  In  the  hyperbolic  method,  which  is  fundamentally 
different  from  conventional  methods,  the  solution  and  the  solution  gradients  are  simultaneously  computed  by 
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solving  a  hyperbolic  system  for  diffusion.  The  hyperbolic  method  was  then  studied  for  advection-diffusion  in 
Ref.  6  with  Residual-Distribution  (RD)  schemes.  Later,  the  method  was  demonstrated  for  the  compressible 
Navier-Stokes  equations  by  a  second-order  Finite- Volume  (FV)  scheme.7  Since  then,  there  have  been  efforts 
to  develop  high-order  hyperbolic  schemes  in  the  finite-volume  method, 8-10  in  the  active  flux  method,11  and 
in  the  RD  method12, 13  for  unsteady  computations. 

In  this  paper,  we  focus  on  the  development  of  hyperbolic  RD  schemes  for  two-dimensional  problems, 
extending  the  previous  work,5,6  with  several  important  advances: 

•  Improved  Accuracy:  We  propose  to  construct  a  second-order  scheme  such  that  it  preserves  exact 
quadratic  solutions,  which  can  be  accomplished  by  the  curvature  correction  technique.14  The  resulting 
scheme  remains  compact,  and  produces  significantly  improved  solution  gradients  over  the  previous 
schemes,  which  do  not  preserve  exact  quadratic  solutions.  Third-Order  Accuracy:  Extending  the  im¬ 
proved  second-order  scheme,  we  construct  a  third-order  scheme  that  preserves  exact  cubic  solutions. 
The  construction  requires  quadratic  least-squares  (LSQ)  gradients  and  a  high-order  source  term  dis¬ 
cretization. 

•  Nonlinear  Equation:  The  improved  schemes  are  extended  to  a  nonlinear  advection-diffusion  equation 
by  the  preconditioned  conservative  formulation  introduced  in  Ref.  7. 

•  High-Order  on  Linear  Elements:  We  demonstrate  that  the  third-order  scheme  does  not  require  curved 
elements  for  curved  boundary  problems;  it  gives  more  accurate  solution  and  gradients  than  the  second- 
order  scheme  on  the  same  linear  grids.  This  is  a  significant  advantage,  because  most  high-order  methods 
require  curved  geometries  to  be  represented  by  high-order  curved  elements  (see  Ref.  3). 

•  Non-Unified  Approach:  Instead  of  the  fully  integrated  approach  of  discretizing  the  hyperbolic  advection- 
diffusion  system  as  in  Ref.  6,  we  discretize  the  advective  and  diffusive  terms  separately.  This  approach 
will  enable  the  extension  to  the  compressible  Navier-Stokes  equations  for  which  the  eigenstructure  of 
the  whole  system  has  not  been  discovered  yet. 

•  Fully  Implicit  Solver:  We  construct  a  fully  implicit  solver  for  both  second-  and  third-order  schemes. 
For  practical  applications,  explicit  iterations  considered  in  Refs.  5  and  6  are  not  efficient  enough,  and 
a  fully  implicit  solver  is  needed.  The  implicit  solver  is  constructed  by  the  exact  residual  Jacobian  of 
the  second-order  scheme.  It  converges  in  at  most  10  iterations  to  reduce  the  residual  by  10  orders  of 
magnitude  for  both  second-  and  third-order  schemes. 

We  demonstrate  these  features  for  a  series  of  test  problems  involving  fully  irregular  isotropic  and  anisotropic 
triangular  grids  and  curved  boundaries.  This  work  serves  as  a  basis  for  the  development  of  high-order  mul¬ 
tidimensional  hyperbolic  RD  schemes  for  more  complex  equations  such  as  the  Navier-Stokes  equations.  The 
extension  of  the  proposed  RD  schemes  to  Navier-Stokes  equations  and  problems  with  shocks  and  disconti¬ 
nuities  will  be  addressed  in  subsequent  papers. 

The  paper  is  organized  as  follows.  In  the  next  section,  we  describe  the  basics  of  the  hyperbolic  RD  scheme: 
residual  evaluations,  boundary  conditions,  and  an  implicit  solver.  In  Section  III,  we  describe  the  construction 
of  the  Low-Diffusion-A  (LDA)  and  the  Streamwise-Petrov-Galerkin  (SUPG)  distribution  matrices  based  on 
a  non-unihed  approach,  in  which  the  advective  and  diffusive  terms  are  treated  independently.  In  Section 
IV,  we  discuss  the  accuracy  issue  of  the  previous  schemes,  and  propose  a  guiding  principle  for  constructing 
improved  schemes.  A  new  second-order  scheme  is  presented  in  Section  V,  followed  by  the  extension  to 
third-order  in  Section  VI.  The  extension  to  nonlinear  advection-diffusion  equations  is  explained  in  Section 
VII.  Numerical  results  are  then  presented  in  Section  VIII  for  both  linear  and  nonlinear  problems,  followed 
by  concluding  remarks  in  Section  IX.  For  completeness,  Fourier  analysis  of  the  constructed  linear  system  is 
also  presented  in  Appendix  A. 

II.  Baseline  RD  Scheme  for  Advection-Diffusion 

This  baseline  scheme  is  the  basis  for  the  development  of  improved  second-order  and  third-order  schemes 
discussed  later.  For  the  purpose  of  our  discussion  and  simplicity,  we  first  describe  the  details  for  the  linear 
advection-diffusion  equation.  Extension  to  the  nonlinear  equation  is  discussed  in  Sec.  VII. 
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A.  Hyperbolic  Advection-Diffusion  System 

Consider  a  two-dimensional  advection-diffusion  equation: 


dtu  +  adxu  +  bdyU  =  v  (dxxu  +  dyyu)  +  s(x,y,u),  (1) 

where  a  and  b  are  constant  advection  speeds,  respectively,  in  x  and  y  directions,  v  (>  0)  is  the  diffusion 
coefficient,  and  s(x,y,u)  is  a  source  term.  Following  Refs.  6  and  12,  we  rewrite  the  above  equation  as  a 
first-order  hyperbolic  advection-diffusion  system: 


dTu  +  adxu  +  bdyU  = 

V  {dxp  +  dvq)  +s(x,y,u), 

(2) 

II 

(dxu  —  p)/Tr, 

(3) 

II 

(dyu  —  q)/Tr, 

(4) 

where  r  is  a  pseudo  time,  Tr  >  0  is  the  relaxation  time,  and  s(x,y,u )  is  a  sum  of  s(x,  y,u)  and  a  physical 
time-derivative  term  discretized  by  an  implicit  time-stepping  scheme  (see  Refs.  12  and  13  for  more  details). 
In  the  vector  form,  the  system  is  written  as 


<9U  .  <9U  dU 

^+A^x+B^y 

where 


u 

a 

-v  0  ’ 

b 

0  -v  ' 

S 

U  = 

P 

,  A  = 

—  1  /Tr 

0 

0 

,  B  = 

0 

0 

0 

,  Q  = 

~P/Tr 

q  _ 

0 

0 

0 

—  1  /Tr 

0 

0 

—q/Tr  _ 

(5) 


(6) 


The  system  is  hyperbolic  in  the  pseudo  time  as  shown  in  Refs.  6.  Towards  the  pseudo  steady  state,  the 
variables  p  and  q  approach  the  solution  gradients  ux  and  uy  and  hence  the  above  equation  becomes  identical 
to  the  original  advection-diffusion  equation  with  the  physical  time-derivative  discretized  by  an  implicit  time¬ 
stepping  scheme  (see  Ref.  12).  This  is  true  for  any  Tr,  and  thus  Tr  is  determined  not  by  physical  constraints 
but  by  optimal  steady  convergence  criteria,6  leading  to  a  non-stiff  hyperbolic  formulation  for  diffusion.  Note 
that  simply  dropping  the  pseudo  time  derivative  will  also  recover  the  original  equation.  In  this  paper,  we 
focus  on  steady  problems,  but  the  resulting  steady  schemes  can  be  made  time  accurate  by  including  a  physical 
time  derivative  term  in  the  source  term,  s(x,y1u),  as  described  in  details  in  Refs.  12  and  13.  In  either  case, 
an  efficient  steady  solver  is  required,  and  its  development  is  one  of  the  objectives  of  the  present  paper. 


B.  RD  Scheme:  Cell-Residual,  Distribution  Matrix,  Nodal  Residual 

We  discretize  the  hyperbolic  advection-diffusion  system  on  unstructured  triangular  grids.  The  domain  is 
divided  into  a  set  {E}  of  arbitrary  triangular  cells  (or  elements),  and  an  associated  set  {J}  of  nodes  (or 
vertices).  The  total  number  of  nodes  is  denoted  by  N.  We  store  the  solutions  ( Uj ,  pj ,  qj )  at  each  node 
j  £  {J}.  In  the  RD  method,  we  first  evaluate  the  residuals  over  the  elements  and  then  distribute  the 
residuals  to  the  nodes  with  a  distribution  matrix.  The  cell-residual  over  an  element  E  £  {E}  is  defined  as 


3 


Figure  1.  Schematic  of  the  residual  distribution  to  local  nodes  and  definition  of  the  inward  unit  normals 
(not-to-scale) . 
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/  (— AUX  —  BUy  +  Q)  dx  dy. 
J  E 


(7) 
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We  assume  a  piecewise  linear  variation  of  U  (i.e. ,  Co  +  cx{x  —  xc )  +  cy(y  —  yc ))  over  the  element,  which 
interpolates  the  three  nodal  solutions,  and  perform  the  integration  to  obtain 

3 

=  -  J2  K*u*  +  QEdnE ,  (8) 

2  =  1 

where  i  is  the  local  vertex  counter  (i.e.,  i  =  1,  2,  3)  for  the  element  E,  d£lE  is  the  element  area,  and 

1  1  3 
Ki  = -(Ahix +Bniy)\ni\  = -Ani\ni\,  QS  =  i^Qi:.  (9) 

i—1 


Note  that  =  (n,;.c ,  nly )  is  the  scaled  inward  normal  to  the  face  (edge)  opposed  to  the  vertex  i,  n,  = 
(see  Fig.  1),  and  n,  =  is  the  unit  inward  normal.  We  note  that  for  every  element 

ELin>  =  o  and  therefore  we  have  K,  =  0. 

The  next  step  is  to  distribute  the  residuals  to  the  three  vertices  by  the  distribution  matrices,  Bf ,  Bf', 
Bf,  as  illustrated  in  Fig.  1.  The  discussion  of  the  distribution  matrix  is  one  of  the  key  contributions  of  the 
present  paper,  but  we  leave  the  choice  open  at  this  point.  We  only  mention  here  that  Bf'  is  a  3x3  matrix, 
which  sums  up  to  the  identity  matrix  over  the  element  for  conservation.  The  distribution  process  results  in 
the  nodal  residual  at  node  j: 


ReSj 


1 

dfl^ 


E&.E 


Bf  $ 


Ee{Ej} 


(10) 


where  dflj  is  the  median  dual  volume  around  the  node  j  (see  Fig.  2),  and  {Ej}  is  a  set  of  triangular  elements 
sharing  the  node  j.  The  nodal  residual  is  an  approximation  of  the  spatial  part  of  the  target  equations,  and 
thus  it  leads  to  a  semi-discrete  form: 

It  can  be  integrated  to  the  steady  state  by  explicit  pseudo-time  stepping  schemes,  as  in  Refs.  5  and  6,  which 
is  significantly  faster  than  conventional  explicit  schemes  because  the  hyperbolic  formulation  eliminates  a 
typical  diffusion  constraint,  0(h2),  where  h  is  a  representative  mesh  spacing,  on  the  explicit  time  step  5  and 
6.  However,  it  still  requires  a  large  number  of  iterations,  especially  for  anisotropic  grids.  To  improve  the 
convergence,  we  develop  an  implicit  solver  in  this  work,  which  is  explained  in  the  next  section. 


Figure  2.  Schematic  of  a  median  dual  volume  around  the  node  j. 


C.  Implicit  Solver 

To  solve  Eq.  (11)  for  the  pseudo  steady  state,  we  drop  the  pseudo-time  derivative,  and  define  the  global 
system  of  steady  residual  equations, 

0  =  Res,  (12) 
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which  consists  of  the  right  hand  side  of  Eq.  (11)  for  all  nodes.  Note  that  the  resulting  residual  equation  has 
become  consistent  with  the  advection-diffusion  equation  (1)  because  the  pseudo-time  derivative  has  been 
dropped.  We  solve  the  system  by  Newton’s  method: 

U'+1  =  U'  +  AU1,  (13) 


where  U  =  (iti,pi,<7i,  U2,P2,q2,  ■  un-,Pn,Qn)  and  l  is  the  iteration  counter.  The  correction  AU!  = 
U 1+1  -  U'  is  determined  as  the  solution  to  the  linear  system: 


<9Res 

<9U 


AU'  =  Res'(U'), 


(14) 


where  Res'  is  the  steady  residual  vector  evaluated  by  U'.  The  Jacobian  matrix  <9Res/dU  is  exact  and 
sparse  because  the  spatial  discretization  is  compact.  For  each  node  j  £  { J},  it  involves  (k  +  1)  3x3  blocks, 
where  k  is  the  number  of  immediate  neighboring  nodes  to  the  node  j.  For  example,  k  =  6  for  the  node  j 
shown  in  Fig.  2  and  therefore,  for  this  particular  node,  there  are  seven  3x3  blocks.  For  a  node  j,  we  have 

k 

-JjAU '  -  £  =  Resj(U'),  (15) 

m—1 


where 


T  dR.es,  T  dR.es, 

~  au,  ’  Jj’m  “  dum  ’  (ibj 

We  may  analytically  evaluate  the  diagonal  and  off-diagonal  entries  of  the  Jacobian  matrix,  i.e. ,  J,  and  Jym> 
but  this  process  is  rather  tedious  for  general  time-dependent  advection-diffusion  problems,  and  even  more 
difficult  for  complex  systems  such  as  the  Navier-Stokes  equations.  To  overcome  the  difficulty,  we  implemented 
an  Automatic  Differentiation  (AD)  tool  based  on  an  operator-overloading  technique  to  evaluate  the  exact 
Jacobians  numerically  through  chain  rule.  Thus,  the  Jacobian  matrix  is  exact  up  to  the  round-off  error 
for  the  baseline  and  second-order  schemes.  The  linear  system  is  relaxed  by  the  sequential  Gauss-Seidel 
relaxation  with  an  under-relaxation  parameter  introduced  for  stabilization  based  on  Fourier  analysis  of  the 
linear  system  (see  Appendix  A  for  more  details  on  stability  analysis)  .  Typically,  the  relaxation  is  performed 
until  the  linear  residual  is  reduced  by  five  orders  of  magnitude  with  a  maximum  relaxation  steps  of  1000. 


D.  Implicit  Boundary  Condition 

Here,  we  discuss  two  Boundary  Condition  (BC)  types:  Dirichlet  (i.e.,  u  is  known),  and  Neumann  (i.e.,  dnu 
is  known).  We  note  that  because  gradients  are  also  the  primary  variables  in  the  hyperbolic  method  (i.e.,  not 
reconstructed  from  the  solution  variable),  the  Neumann  type  BC  is  treated  similarly  as  the  Dirichlet  type 
BC. 


Figure  3.  Schematic  of  boundary  nodes  for  formulation  of  implicit  boundary  condition. 


For  a  Dirichlet  type  BC,  the  following  set  of  equations  is  imposed  at  the  boundary  nodes: 


u  —  Ub 

(p,  q)  ■  h 

(■ dxu  -  p,  dyu  -  q)  -  fib 


0, 

0, 


(17) 

(18) 
(19) 


where  Ub  is  the  given  boundary  value  and  dsUb  is  the  derivative  of  the  Ub  along  the  boundary,  which  can 
be  computed  with  the  given  Ub ■  The  fib  and  tb  are,  respectively,  the  unit  normal  and  tangent  vectors  at 
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the  boundary  nodes  (see  Fig.  3).  For  a  strong  Dirichlet  type  BC,  we  specify  the  BC  value  at  the  boundary 
nodes  and  solve  the  hyperbolic  system  accordingly.  Here,  we  discretize  and  solve  Eq.  (19)  according  to 
the  hyperbolic  RD  scheme  described  in  the  previous  section.  That  is,  the  residual  vector  at  the  node  jb  is 
redefined  as 


'U'jb 


(. Pjb  ’  Qjb  )  *  ds'U'b 


(20) 


L  ResJb(2)«6x  +Resib(3)«by  J 

where  ReSj6(2)  and  ResJ6(3)  are  the  nodal  residuals  of  the  second  and  third  equations  of  the  first-order 
hyperbolic  system  (i.e. ,  Eq.  3  and  Eq.  4,  respectively),  computed  as  described  in  the  previous  sections. 

For  a  Neumann  type  BC,  we  impose  the  following  equations  to  the  boundary  nodes: 


dTu  +  adxu  +  bdyU 

=  v(dxp  +  dyq)  +  s, 

(21) 

( P >  q)  ■  nb 

=  dnub, 

(22) 

)xu  -  p,  dyu  -  q)  -tb 

=  o, 

(23) 

where  dnUb  is  the  given  gradient  of  the  u  on  the  boundary  in  the  normal  direction,  and  therefore  is  known. 
For  the  strong  formulation  of  a  Neumann  type  BC,  we  discretize  and  solve  Eqs.  (21)  and  (23)  according  to 
the  RD  scheme;  that  is 


ResJ6(l) 


( Pjh,qjb )  ’  nb  -  dnub 


(24) 


ResJb(2)ibi:  +Hesjb(3)tby 


where  ResJb(l)  is  the  nodal  residual  of  the  first  equation  of  the  first-order  hyperbolic  system  computed  as 
described  in  the  previous  sections. 

In  either  of  the  Dirichlet  or  the  Neumann  type  BC,  the  modified  residual  is  incorporated  into  the  implicit 
solver  with  the  exact  Jacobians.  In  doing  so,  a  care  must  be  taken  to  avoid  zero  diagonal  entries  in  the 
Jacobian;  the  second  and  third  components  in  the  boundary  residual  need  to  be  exchanged,  depending  on 
the  magnitude  of  hbx  and  hb  ,  to  avoid  zero  diagonals. 


E.  Remark  on  variable  advection  coefficients 


If  the  advection  coefficients  a  and  b  are  functions  of  (x,y),  advective  fluxes  are  not  necessarily  conservative. 
In  this  case,  we  estimate  the  advection  vector  by  the  arithmetic  average  over  an  element  to  write 


Ih  + 


=dU  ^ 
+  B— —  —  Q> 
oy 


(25) 


where  A  and  B  are  defined  with  the  averaged  advection  vector  (a,  b): 


a 

-v  0  " 

b 

0  -v  ] 

A  = 

-1/Tr 

0 

0 

,  B  = 

0 

0 

0 

0 

0 

0 

.  -1/Tr 

0 

0 

The  baseline  RD  scheme  is  then  applicable.  If  conservative  fluxes  exist,  then  the  nonlinear  formulation 
described  in  Sec.  VII  can  be  employed  to  construct  an  RD  scheme. 


III.  Distribution  Matrices:  Non-Unified  Approach 

The  RD  scheme  is  defined  by  the  combination  of  the  cell-residuals  and  the  distribution  matrix.  The 
distribution  matrix  is  mainly  responsible  for  the  dissipative  behavior  and  the  stability  while  the  cell-residuals 
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determine  the  order  of  accuracy.  Unbounded  distribution  matrices  can  also  affect  the  order  of  accuracy,  but 
we  consider  only  bounded  ones  in  this  paper.  The  construction  of  the  distribution  matrix  often  requires 
the  knowledge  of  the  eigenstructure  of  the  target  system.  The  eigenstructure  of  the  hyperbolic  advection- 
diffusion  system  is  known6  and  various  distribution  matrices  can  be  constructed.  However,  the  extension 
to  the  compressible  Navier-Stokes  equations  is  not  possible  at  the  moment  because  the  eigenstructure  has 
not  been  successfully  worked  out  yet.  In  order  to  develop  hyperbolic  RD  schemes  that  can  be  extended  to 
the  compressible  Navier-Stokes  equations,  we  consider  a  non-unified  approach  introduced  for  FV  schemes 
in  Refs.  7  and  9,  respectively,  for  the  compressible  Navier-Stokes  and  advection-diffusion  equations.  In  this 
approach,  we  construct  distribution  matrices  for  advection  and  diffusion  separately  and  then  combine  them 
to  derive  a  distribution  scheme  for  the  advection-diffusion  equation.  This  approach  easily  extends  to  the 
compressible  Navier-Stokes  equations  because  the  eigenstructures  of  the  inviscid  terms  and  the  hyperbolized 
viscous  terms  can  be  fully  analyzed  independently.7 


A.  Unified  Approach 

To  illustrate  the  point  of  the  non-unified  approach,  we  first  recall  the  unified  approach  of  Ref.  6.  Consider 
an  arbitrary  unit  vector  h  =  (nx.  ny ) .  The  unified  advection-diffusion  flux  Jacobian  is  defined  as: 


A  n  =  A  nx  +  B  fiy  = 


an  -vnx  -vny 
— nx/Tr  0  0 

—riy  /Tr  0  0 


(27) 


where  an  =  ahx  +  bhy  is  the  advection  velocity  projected  onto  the  unit  vector  n.  The  Jacobian  An  has  the 
following  real  eigenvalues: 


Ai  = 


1 

/  „ 

,  1 

/  „  4i/  ' 

2 

dn  * 

V  "  +  Tr  _ 

’  Aa~2 

/a”  +  Tr  _ 

A;!  —  0, 


where  Ai  <  0  and  A2  >  0.  The  arbitrary  but  positive  relaxation  time,  Tr ,  is  defined  as 

_  ... 

T  = 

r  — 


V a2  +  b2  +  v/Lr  ’ 


(28) 


(29) 


where  the  length  scale  Lr  is  a  quantity  of  0(1),  which  may  be  taken  as  Lr  =  l/2ir,  as  used  here,  or  from  an 
optimal  formula  as  derived  in  Ref.  6.  The  right  and  left  eigenvectors  for  the  Jacobian  An  are  given  by 


(30) 


1 

> 

■i 

— A2  Tr 

0 

-1  /Tr 

—  A2  fly 

R„  = 

nx 

nx 

-fry 

T  1 

’  ~  A1-A2 

1  /Tr 

^inx 

Al  fly 

fl.y 

fly 

nx 

0 

—  (Ai  —  A2 )ny  (Ai  —  \2)nx  J 

The  unified  Jacobian  An  can  be  rewritten  based  on  its  right  and  left  eigenvectors  as 


-A-rc  BnAnLn, 


(31) 


where  A  is  the  diagonal  matrix  with  the  entries  of  the  eigenvalues  defined  in  Eq.  (28).  The  Jacobian  matrix 
can  be  decomposed  as  6 


where 


A1II1  +  A^ 

>n2  —  An 

+  A  +  , 

(32) 

Ai 

Ai  A2  Trfax 

Ai  A2  TrTly 

fax /Tr 

-\2fal 

X2  fax  fay 

,  (33) 

/Tr 

^2  fax  fay 

-A  2nl 

~  A2 

-  Ai  A2  Trfax 

—  \l\2Trfay 

1 

fix  /Tr 

\ifa2x 

X\faXfay 

(34) 

fay  /Tr 

\\faXfay 

A  iK  . 
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Thus,  the  K;  defined  in  Eq.  (9)  can  be  written  as 


k,  =  k; 


■K+ 


(35) 


where 


K- 

K+ 


2Ai,i  Hi,*  |n,|  —  2^n 


n,; 


2^2,i  n2)i  |n» 


1 


=  2  A™»  I11*!' 


(36) 

(37) 


These  matrices  are  often  required  for  the  construction  of  the  distribution  matrix,  and  as  can  be  seen  from 
the  above,  they  require  the  eigenvalues  and  the  eigenvectors,  which  are  not  known  at  present  for  a  hyperbolic 
formulation  of  the  compressible  Navier-Stokes  equations. 


B.  Non-Unified  Approach 

In  the  non-unihed  approach,7,9  we  treat  the  advective  and  diffusive  terms  separately.  Therefore  we  rewrite 
Eq.  (5)  as 

<9U  9F  dG  „ 

^+dx+^y~Q’  (38) 

where 


p  _  pQ 


.  Fd  = 


G  =  Ga  +  G“  = 


au 

—up 

0 

+ 

—u/Tr 

0 

0 

bu 

-uq 

0 

+ 

0 

0 

1 

s 

1 

_ i 

(39) 


(40) 


Similarly,  the  flux  Jacobian  Ara  is  decomposed  into  two  advective  and  diffusive  fluxes,  A“  and  A*,  respec¬ 
tively: 

A  —  A°  4-  \d 


where 


Q"n 

0 

0  ’ 

0 

—ufix 

—  Ufly 

A  a  — 

0 

0 

0 

A  d  — 

5 

fix /Tr 

0 

0 

0 

0 

0 

1 

3> 

0 

0 

(41) 


(42) 


The  distribution  matrix  is  constructed  separately  for  the  advective  and  diffusive  terms  based  on  the  corre¬ 
sponding  flux  Jacobian,  and  then  the  two  matrices  are  combined  to  form  a  matrix  for  the  advection-diffusion 
equation.  This  approach  can  be  extended  to  the  compressible  Navier-Stokes  equations  because  the  eigen- 
structure  of  the  hyperbolized  viscous  terms  is  available. ' 

The  advective  Jacobian  has  only  one  non-zero  eigenvalue  (Aa  =  an)  and  thus  we  can  simply  decompose 
it  according  to  the  sign  of  Aa: 


4“  =  ia  +  A“  -=■ 


max(0,  an ) 

0 

0  " 

min(a„,0) 

0 

0  ’ 

0 

0 

0 

+ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(43) 


The  diffusive  Jacobian  has  two  non-zero  eigenvalues,  Af  =  —  yJu/Tr,  Af[  =  +yJv/Tr,  with  the  right  and  left 
eigenvectors  formally  in  the  form  of  Eq.  (30).  Note,  however,  that  they  are  numerically  different  because  the 
eigenvalues  are  different.  The  diffusive  Jacobian  can  be  split  as 


K  =  «K  =  A?m  +  A*n2  =  a: 


d+ 


(44) 
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where  FEi  and  n2  are  given  by  Eq.  (33)  and  Eq.  (34),  respectively,  with  Ai  =  Af  and  A2  =  Aff. 
We  define  Kf  by 


K?=Ul|ni|=Kf +K 


where  the  negative  and  positive  matrices,  Kf  and  Kf+  are  given  by 


d+ 


(45) 


Kd 

= 

=  Iwl, 

(46) 

Kd+ 

2 

=  n' 

II 

tO  1  H 

> 

p 

(47) 

Various  distribution  matrices  can  be  employed  for  the  hyperbolic  RD  schemes:  Lax-Wendroff  (LW),5 
Low-Diffusion- A  (LDA),6  and  Streamline-Upwind-Petrov-Galerkin  (SUPG)  distribution  matrices.  In  the 
following  sections,  we  describe  the  LDA  and  the  SUPG  distribution  schemes  and  their  formulations  in  the 
context  of  hyperbolic  RD  schemes  (the  LW  formulation  is  similar  to  the  SUPG).  For  completeness,  both 
unified  and  non-unified  approaches  are  presented. 


C.  Distribution  Matrix:  Low-Diffusion-A  (LDA) 

The  LDA  scheme  is  a  multidimensional  upwind  scheme,  meaning  that  B^13*  becomes  zero  when  no  signal 
is  being  sent  to  the  the  upwind  nodes.  In  the  unified  approach,  the  distribution  matrix  B  f  for  the  LDA 
scheme15, 16  is  defined  as 


B)DA  =  K 


(48) 


For  the  non-unified  approach,  the  distribution  matrix  should  separately  be  constructed  with  the  K)5  °  and 
the  K)1" d  terms,  corresponding  to  the  advection  and  the  diffusion  components,  respectively.  The  diffusion 
component  of  the  distribution  matrix,  Bf,  is  similar  to  Eq.  (51),  except  that  the  K)1  is  replaced  with  the 

K+d: 

-l 


B“  =  K 


+  d 


(49) 


Note  that  the  inverse  matrix  is  usually  evaluated  numerically,  as  done  here,  but  the  analytical  expression  for 
Bf  is  also  available  in  Ref.  17,  which  we  used  it  only  for  verification  purposes.  The  distribution  matrix  for  the 
advection  term,  B“,  reduces,  however,  to  a  scalar  equation  as  given  below  because  there  is  no  contribution 
from  the  gradient  equations: 


B“ 


/?2lda  0  0 

0  0  0 
0  0  0 


ftDA 


max(0,  ani)\rii\ 
j=i  m&x(0,  anj )  |  n.j  \ 


(50) 


These  two  matrices  are  combined  in  the  following  form  to  define  the  LDA  scheme  for  the  advection-diffusion 
equation: 


B) 


=  (I  -  IW)B“  +  IWB i 


(51) 


where 


to  0  0 

0  10, 

0  0  1 


(52) 


and  w  is  a  weighting  function.  From  a  detailed  analysis  of  a  one-dimensional  hyperbolic  advection-diffusion 
system,  the  weighting  function  can  be  evaluated  as6 


to  = 


2 

Re +  2’ 


(53) 
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where  we  define  the  Re  for  the  two-dimensional  system  as 


Re  =  \J  0?  +  b2/v.  (54) 

The  weights  have  been  introduced  for  two  purposes.  One  is  to  guarantee  the  conservation: 

EB"DA  =  i’  (55) 

and  the  other  is  to  ensure  that  the  scheme  approaches  a  pure  advection  scheme  in  the  limit  of  Re  — >  oo, 
and  a  pure  diffusion  scheme  as  Re  0.  Note,  for  example,  that  /3fDA  is  0(1)  even  for  a  vanishingly  small 
advective  vector,  and  thus  it  will  affect  the  distribution  even  in  the  diffusion  limit  if  the  weighing  function 
is  independent  of  Re. 

In  our  experience,  hyperbolic  RD  schemes  with  the  LDA  distribution  matrix,  constructed  by  the  non- 
unified  approach,  have  an  instability  problem  in  the  linear  relaxation  on  some  anisotropic  triangular  grids 
even  with  the  employment  of  uder-relaxation  parameter  obtained  from  Fourier  analysis  (see  Appendix  A). 
On  the  other  hand,  hyperbolic  RD  schemes  can  be  successfully  stabilized  with  a  suitable  under-relaxation 
parameter  in  the  case  of  the  SUPG  distribution  matrix,  which  is  discussed  next. 


D.  Distribution  Matrix:  Streamlme-Upwind-Petrov-Galerkin  (SUPG) 

The  finite-element  SUPG  scheme18  is  known  to  be  applicable  to  the  RD  framework.19  In  the  unified  approach, 
the  SUPG  distribution  matrix  for  triangular  grids  is  given  by 


B?UPG=  3I- 


K, 


w=i 


(56) 


which  consists  of  the  Galerkin  part  (the  first  term)  and  the  stabilization  part  (the  second  term).  In  the 
non-unified  approach,  we  construct  the  SUPG  distribution  matrix  as  follows: 

Bf PG  =  *1  +  D“  +  Df ,  (57) 

where  D“  and  Df  are  the  stabilization  terms  defined  independently  for  the  advective  and  diffusive  terms, 


D 


dfPG  0  0 
0  0  0 
0  0  0 


JSUPG 


1  _ _ 

2  X^f=i  max(0,  ani)\ni\ 


D“  = 


(58) 


(59) 


Note  that  the  denominator  of  Eq.  (58)  cannot  vanish  unless  the  advection  vector  is  exactly  zero.  A  small 
numerical  value  of  the  order  of  machine  zero  may  be  added  to  the  denominator  in  order  to  avoid  zero  division 
in  the  case  of  zero  advection  vector.  Another  possibility  is  to  completely  remove  D';‘  term  for  zero  advection 
vector;  both  approaches  yield  identical  results. 

It  is  also  important  to  note  that  we  have  Bf  =  JU=1  Bf  =  0,  and  therefore  the  conservation 
property  is  guaranteed: 

EB?UPG  =  L  (6°) 

i= 1 

A  weighting  function  can  be  introduced  in  the  stabilization  terms  of  Eq.  (57), 

BfUPG  =  ^I  +  (l  —  w)D?+wD?,  (61) 

where  oj  is  defined  such  that  w  — >  0  as  Re  — >  oo ,  and  u  — >  1  as  Re  — >  0,  so  that  the  scheme  will  properly  reduce 
to  pure  advection  and  diffusion  schemes  as  Re  — »  oo  and  Re  — >  0,  respectively.  The  weighting  function  oj 
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can  be  evaluated,  based  on  a  detailed  analysis  for  a  one-dimensional  hyperbolic  advection-diffusion  system,6 
as 


u  = 


2 

i?e  +  2  ’ 


(62) 


where  we  define  the  Re  for  the  two-dimensional  system  as 


Re  =  \J  a2  +  b2 / v. 


(63) 


In  our  study,  however,  we  observed  no  noticeable  change  in  the  results  between  Eq.  (57)  and  Eq.  (61)  with  the 
above  weighting  function.  Therefore,  the  results  shown  in  this  paper  are  all  based  on  the  SUPG  distribution 
matrix  without  the  weighting  function.  Details  of  the  effects  of  the  weighting  function  are  a  subject  for 
future  study.  We  remark  that  unlike  in  the  LDA  formulation,  where  the  weighting  function  is  required  to 
preserve  conservation,  the  SUPG  formulation  does  not  require  any  weighing  function  because  conservation 
is  guaranteed  regardless. 

As  can  be  seen  from  above,  the  SUPG  distribution  matrix  is  constructed  based  on  a  decomposition  of 
the  stabilization  term  into  advective  and  diffusive  stabilization  terms.  Therefore,  the  SUPG  distribution 
matrix  can  be  extended  to  a  hyperbolic  formulation  of  the  compressible  Navier-Stokes  equations,  for  which 
the  inviscid  and  viscous  terms  can  be  analyzed  independently.7 


IV.  New  Design  Principle  for  Advection-Diffusion  Equation 

An  important  feature  of  the  RD  schemes  critical  to  the  accuracy  condition  is  the  ability  to  preserve 
exact  polynomial  solutions  on  arbitrary  grids.  Historically,  RD  schemes  have  been  designed  to  preserve 
linear  exact  solutions  for  hyperbolic  systems  of  conservation  laws.  This  is  called  the  linearity  preservation. 
Such  schemes  are  known  to  yield  second-order  discretization  errors.  The  baseline  hyperbolic  RD  scheme5 
is  constructed  based  on  the  linearity  preserving  property.  Although  it  gives  second-order  accuracy  for  all 
variables  on  regular  triangular  grids,5  it  was  found  later  that  the  order  of  accuracy  of  the  gradient  variables 
(i.e.,  p  and  q)  deteriorates  to  first-order  on  irregular  grids.  In  a  subsequent  work  presented  in  Ref.  6,  a  better 
scheme  was  proposed  that  attempted  to  improve  the  accuracy  of  the  gradients  by  upgrading  the  evaluation 
of  the  second  and  third  components  of  the  residual.  In  this  paper,  the  scheme  of  Ref.  6  is  referred  to  as  the 
RD-JCP2010  scheme.  Results  given  in  Ref.  6  showed  that  the  RD-JCP2010  scheme,  relative  to  the  baseline 
scheme,  improves  the  order  of  accuracy  of  the  gradient  variables,  but  no  results  were  presented  for  the  quality 
of  the  predicted  gradients.  As  will  be  shown  later  in  this  paper,  we  found  that  for  a  different  test  problem 
and  a  more  general  exact  solution,  the  same  RD-JCP2010  scheme  does  not  yield  second-order  accuracy  for 
the  gradient  variables  (i.e.,  p  and  q),  and  also  generates  severe  oscillations  in  these  variables  within  the 
domain.  It  appears  that  the  discretization  of  the  source  terms  (including  those  arising  from  the  hyperbolic 
formulation)  was  not  fully  compatible  with  the  flux  terms,  and  it  is  not  clear  to  us,  at  the  moment,  how  to 
improve  the  source  term  discretization  to  develop  a  fully  second-order  scheme.  To  overcome  this  difficult 
problem,  we  instead  propose  to  improve  the  flux  terms  such  that  the  residual  vanishes  for  quadratic  solutions 
for  the  advection-diffusion  equation.  The  exactness,  as  will  be  demonstrated  later,  has  a  significant  impact 
on  the  accuracy  and  the  quality  of  the  gradient  variables  on  irregular  grids.  The  construction  of  such  schemes 
is  very  difficult  in  general  because  second-order  gradients  need  to  be  recovered  from  second-order  accurate 
solutions  on  irregular  grids  for  the  diffusive  term.  However,  the  construction  is  quite  straightforward  in  the 
hyperbolic  method  since  the  gradients  are  computed  simultaneously  with  the  solution  variable  u.  The  scheme 
also  extends  easily  to  the  third-order  by  imposing  the  exactness  for  cubic  solutions.  Note  that  the  exactness 
is  not  a  necessary  condition  for  accuracy.  For  example,  the  third-order  RD  scheme  proposed  in  Ref.  20 
does  not  preserve  cubic  solutions,  but  was  shown  to  produce  third-order  accuracy  on  smooth  unstructured 
triangular  grids. 

In  the  next  two  sections,  we  describe  how  to  construct  second-  and  third-order  accurate  schemes  that, 
respectively,  preserve  quadratic  and  cubic  solutions  on  arbitrary  triangular  grids. 
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V.  Improved  Second-Order  Scheme:  RD-CC2 


Consider  the  cell-residuals  for  the  baseline  hyperbolic  RD  scheme: 


®u=J  (-aux  -  buy  +  upx  +  uqv  +  s)dxdy  =  ^  ^  {-aniUi\n,\  +  i/pjnix  +  vqiniy )  +  sEdflE 


=  jr  I  {ux~  p)dxdy  =  |  )  . 
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Tr  \  2 


»= l 


=  - 
9  Tr 


[  ( uy  -  q)dxdy  =  A  (  *  V 

lr  \ 2  i= 1 


(64) 

(65) 

(66) 


These  residuals  are  constructed  based  on  the  assumption  that  all  variables  (i.e.,  u,  p,  q)  are  linear  over  the 
element,  and  therefore  they  do  not  preserve  quadratic  solutions.  Note  that  quadratic  solutions  imply  linear 
gradients  (i.e.,  ux,  uy,  p,  q),  and  the  above  cell  residuals  are  already  exact  for  linear  gradients.  Therefore, 
to  preserve  quadratic  solutions,  the  integrals  of  ux  and  uy  terms  need  to  be  exact  for  quadratic  solutions. 
A  simple  way  to  accomplish  this  is  to  reconstruct  a  quadratic  element  by  interpolating  the  solution  at  the 
midpoint  of  each  side  in  the  element  (see  Fig.  4).  The  midpoint  value  umi  is  estimated  by  the  Hermite 


3 


Figure  4.  Reconstructed  P2  element  with  virtual  midpoints. 


interpolation  along  the  side: 


umi  =  Ui  -  -  (ApiAxi  +  AqiAyi) , 

O 


(67) 


where  Tii  is  the  arithmetic  average  of  the  solution  values  at  the  two  end  nodes,  and  A( ),  denotes  the  difference 
of  the  nodal  values  taken  counterclockwise  along  the  edge  opposite  to  the  node  i,  e.g.,  Ax 3  =  x\  —  X2-  Note 
that  the  midpoint  value  umi  is  exact  for  quadratic  solutions  (and  linear  gradients).  Once  the  quadratic 
element  is  reconstructed,  the  cell-residuals  are  evaluated  as  line  integrals  with  Simpson’s  rule  applied  along 
each  side.  This  reconstruction  approach  was  first  introduced  in  Refs.  21  and  22,  and  later  employed  in 
the  form  of  the  Green-Gauss  gradient  with  a  high-order  curvature  correction  term.6, 14,23  In  this  work,  we 
employ  the  latter  implementation  to  upgrade  the  baseline  residuals.  Note  that  the  formula  in  Eq.  (67) 
usually  requires  LSQ  gradient  reconstruction  to  obtain  the  nodal  gradients,14, 21-23  but  it  is  not  necessary  in 
the  hyperbolic  method  because  p(=  ux )  and  q(=  uy)  are  already  available  at  nodes  as  a  part  of  the  primary 
variables. 

Applying  the  curvature  correction,  we  modify  the  cell-residuals  of  the  baseline  RD  scheme  as 
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(71) 


where  the  curvature  correction  term,  Sf ,  is  defined  as 

Si  =  \  ( A  pi  Ax.i  +  AqiAyi). 

b 

Note  that  in  Ref.  6,  5“  was  mistakenly  shown  with  a  negative  sign.  One  may  find  it  convenient  that  the 
compact  RD-CC2  scheme  can  be  implemented  as  an  add-on  to  the  baseline  hyperbolic  RD  scheme: 
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The  modified  cell-residuals  are  distributed  by  the  SUPG  distribution  matrix.  The  RD  scheme  based  on 
these  cell-residuals  is  referred  to  as  the  RD-CC2  scheme.  This  scheme  is  compact  because  no  LSQ  gradient 
reconstruction  is  necessary.  The  exactness  for  quadratic  solutions  has  been  confirmed  numerically,  and  are 
compared  with  the  baseline  RD  and  RD-JCP2010  schemes.  These  results  are  summarized  in  Table  1, 


Table  1.  Exactness  of  the  hyperbolic  RD  equations  for  exact  quadratic  polynomial  solutions  on  irregular  grids. 
The  last  column  indicates  whether  the  exactness  (/)  is  based  on  cell  and/or  nodal  residuals. 


Hyperbolic  RD  Scheme 

Quadratic  solution 

u  eq.  p  eq.  q  eq. 

Cell/Nodal 

Baseline 

X 

X 

X 

none 

RD-JCP2010 

X 

/ 

/ 

cell 

RD-CC2 

/ 

/ 

/ 

both 

The  truncation  error  orders  for  the  RD-CC2  scheme,  the  baseline  RD  scheme,  and  the  RD-JCP2010 
scheme  are  provided  in  Table  2.  The  truncation  error  orders  have  been  determined  by  substituting  a  smooth 
exact  solution  (taken  from  Ref.  24)  into  the  residuals  for  a  series  of  irregular  triangular  grids.  Clearly, 
all  schemes  have  the  same  truncation  error  order  for  the  nodal  residual.  A  known  theory  for  hyperbolic 
conservation  laws25,26  states  that  the  discretization  error  of  0(h 2)  leads  to  the  truncation  error  of  0(h). 
However,  our  numerical  results  indicate  that  the  converse  is  not  necessarily  true.  As  we  will  show  later,  only 
the  RD-CC2  scheme,  which  has  0(h)  truncation  order,  gives  a  discretization  error  of  0(h2). 

Table  2.  Comparison  between  numerical  truncation  error  of  the  nodal  residuals  (i.e.,  T.E.  =  JT  Bb#b) 
obtained  with  RD-CC2  and  non-genuine  second-order  RD  schemes  on  irregular  grids. 


Hyperbolic  RD  Scheme 

Nodal  Residuals 

u  equation  p  equation  q  equation 

Baseline 

0(h) 

0(h) 

0(h) 

RD-JCP2010 

0(h) 

0(h) 

0(h) 

RD-CC2 

0(h) 

0(h) 

0(h) 

We  remark  that  the  RD-CC2  scheme  is  different  from  the  RD-JCP2010  scheme,6  which  has  the  curvature 
correction  only  in  and  residuals.  Thus,  the  whole  residual  vector  of  RD-JCP2010  is  not  exact  for 
quadratic  solutions  as  confirmed  in  Table  1.  The  RD-CC2  scheme,  on  the  other  hand,  has  the  curvature 
correction  applied  to  all  the  residuals,  and  therefore  is  exact  for  quadratic  solutions  (and  linear  gradients). 

To  solve  the  discrete  equations,  we  employ  the  implicit  solver  described  in  Sec.  C.  The  RD-CC2  scheme 
is  compact  and  therefore  we  linearize  the  residual  vector  exactly  by  using  the  AD  tool.  The  implicit  solver 
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is  thus  Newton’s  method  and  converges  to  machine  zero  (i.e. ,  10  orders  of  magnitude  reduction  in  all  the 
residuals)  typically  within  3-5  Newton  iterations. 


VI.  Third-order  Scheme:  RD-CC3 


To  extend  the  compact  RD-CC2  scheme  to  a  third-order  scheme,  we  modify  the  cell-residuals  so  that  they 
preserve  exact  cubic  solutions  and  quadratic  gradients.  Following  the  procedure  explained  in  Sec.  V  for  the 
RD-CC2,  and  thanks  to  the  fact  that  both  the  Hermite  interpolation  and  Simpson’s  rule  are  exact  for  cubic 
functions,  we  employ,  again,  the  curvature  correction  technique.  Extension  to  the  third-order,  therefore, 
requires  two  additional  steps:  1)  high-order  source  term  discretization,  and  2)  quadratic  LSQ  gradients  for 
p,  q ,  and  the  source  term  s.  The  high-order  discretization  of  the  source  term  is  a  critical  piece  for  achieving 
third-order  accuracy.  We  derive  a  formula  by  estimating  the  midpoint  value  by  the  Hermite  interpolation 
and  then  exactly  integrating  the  source  terms  assuming  a  quadratic  variation  over  the  element: 


where 


IE 


s  dxdy 
p  dxdy 
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Using  these  high-order  source  discretizations,  we  construct  the  cell-residuals  that  preserve  cubic  solutions 
(and  quadratic  gradients)  as  follows: 
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which  can  also  be  written  in  the  form  of  an  extension  to  the  compact  RD-CC2  scheme: 
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The  curvature  correction  terms  require  quadratic  LSQ  gradients  for  p.  q,  and  s,  but  not  for  u  because 
p(=  Ux)  and  q{=  uy)  are  already  available  at  nodes.  Again,  the  cell-residuals  are  distributed  by  the  SUPG 
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distribution  matrix.  The  hyperbolic  RD  scheme  based  on  these  cell-residuals  is  referred  to  as  the  RD-CC3 
scheme.  The  exactness  for  cubic  solutions  has  been  confirmed  numerically,  as  shown  in  Table  3.  The  results 
of  a  numerical  truncation  error  analysis  for  the  RD-CC3  are  provided  in  Table  4  and  compared  with  the 
truncation  error  of  RD-CC2. 

Table  3.  Exactness  of  the  hyperbolic  RD  equations  for  exact  quadratic  and  cubic  polynomial  solutions  on 
irregular  grids.  The  last  column  indicates  whether  the  exactness  (/)  is  based  on  cell  and/or  nodal  residuals. 


rr  i  i-  t>t,  n  i  Quadratic  solution  Cubic  solution  „.  ,  , 

Hyperbolic  RD  Scheme  Eiemental/Nodai 

u  eq.  p  eq.  q  eq.  u  eq.  p  eq.  q  eq. 

RD-CC2  ///XXX  both 

RD-CC3  //////  both 


Table  4.  Comparison  between  numerical  truncation  error  of  the  nodal  residuals  (i.e.,  T.E.  =  JT  Bb#b) 
obtained  with  the  RD-CC2  and  RD-CC3  schemes  on  irregular  grids. 


Hyperbolic  RD  Scheme 

RD-CC2 

RD-CC3 


Nodal  Residuals 

u  equation  p  equation  q  equation 

0(h)  0{h)  0{h) 

0(h2)  0{h2)  0(h2) 


The  discrete  equations  are  solved  by  the  implicit  solver  described  in  Sec.  C.  We  do  not  attempt  to 
linearize  the  RD-CC3  scheme  exactly,  but  simply  use  the  exact  Jacobian  of  the  compact  RD-CC2  scheme. 
The  implicit  solver  is  therefore  not  precisely  Newton’s  method,  but  it  converges  so  rapidly  that  fully  converged 
numerical  solutions  (i.e.,  10  orders  of  magnitude  reductions  for  all  residuals)  are  obtained  typically  within 
10-15  Newton  iterations  as  will  be  demonstrated  later  in  Sec.  VIII. 

Remark:  The  RD-CC2  and  RD-CC3  schemes  with  the  SUPG  distribution  matrix  can  be  considered  as 
economical  and  powerful  alternatives  to  high-order  finite-element  methods.2'  They  are  economical  mainly  for 
three  reasons:  1)  the  number  of  linear  relaxations  in  the  implicit  solver  increasers  linearly,  not  quadratically 
as  is  typical  for  diffusion  problems,  with  grid  refinement.  2)  the  second-order  advection-diffusion  scheme  is 
compact,  and  thus  allows  an  efficient  construction  of  Newton’s  method.  Furthermore,  the  same  compact 
Jacobian  still  yields  rapid  convergence  (comparable  to  Newton’s  method)  for  the  third-order  scheme.  3) 
second-  and  third-order  accuracy  can  be  achieved  (as  will  be  demonstrated  later  in  Sec.  VIII)  on  linear 
triangular  elements  without  introducing  curved  elements.  Thus,  higlr-order  curved  grids  are  not  needed,  and 
these  schemes  are  directly  applicable  to  existing  grids  composed  of  linear  triangular  elements.  These  schemes 
are  also  powerful  in  that  they  are  capable  of  producing  highly  accurate  and  smooth  gradients,  to  the  same 
order  of  accuracy  as  that  of  the  main  solution  variable  u,  on  isotropic  and  anisotropic  irregular  grids.  These 
advantages  will  be  demonstrated  in  Sec.  VIII. 

VII.  Nonlinear  Advection-Diffusion  Equation 

In  this  section,  we  describe  the  discretization  of  a  nonlinear  hyperbolic  advection-diffusion  equation.  As 
an  example,  we  discuss  the  construction  of  the  RD-CC3  scheme  with  the  SUPG  distribution  matrix.  Note 
that  the  compact  RD-CC2  is  a  subset  of  RD-CC3  scheme  and  therefore  the  discussion  here  also  applies  to 
the  RD-CC2  scheme.  A  similar  procedure  can  be  applied  to  other  distribution  schemes.  Throughout  the 
discussion,  we  only  consider  the  non-unified  approach,  which  is  applicable  to  more  complex  systems  such  as 
the  compressible  Navier-Stokes  equations. 

A.  Nonlinear  Hyperbolic  Advection-Diffusion  Equation 

Consider  the  following  two-dimensional  nonlinear  advection-diffusion  equation: 

dtu  +  dxf  +  dyg  =  dx(vdxu)  +  dy(vdyu )  +  s(x ,  y ,  u),  (87) 
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where  /  and  g  are  nonlinear  functions  of  u,  and  v  =  v(u).  The  advection  speeds  in  x  and  y  directions  are 
therefore  a(u)  =  df  /du  and  b(u)  =  dg/du ,  respectively.  Using  the  preconditioned  formulation  proposed  for 
nonlinear  systems  in  Ref.  7,  we  construct  a  hyperbolic  system  as 


dTu  +  dxf  +  dyg 

=  dxp  +  dyq  +  s, 

(88) 

Tr  o 

— oTp 

V 

=  dxu  -  p/v, 

(89) 

T 

±r  0 

—  Orq 

V 

=  dyu  —  q/ v, 

(90) 

where  the  variables  p  and  q  are,  in  the  pseudo  steady  state,  equivalent  to  the  diffusive  fluxes  in  x  and  y 
directions,  respectively.  As  before,  the  physical  time  derivative  can  be  incorporated  into  the  source  term 
s  as  is  done  in  Ref.  13,  but  we  focus  on  the  steady  equation  here.  Note  that  the  system  reduces  to  the 
target  equation,  i.e. ,  Eq.  (87),  in  the  pseudo  steady  state.  The  system  is  written  in  the  vector  form  as  a 
preconditioned  conservative  equation,  with  the  preconditioning  matrix  P: 
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The  flux  Jacobians  can  be  obtained  by  multiplying  both  sides  of  Eq.  (91)  by  the  preconditioned  matrix 
P,  and  arrive  at 


<9F 
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Note  that  these  preconditioned  Jacobian  matrices  (i.e.,  A  and  B)  are  slightly  different  from  the  ones  obtained 
for  the  linear  hyperbolic  advection-diffusion  system.  Hence,  their  eigen-structures  are  also  different  from 
those  given  for  the  linear  system  in  Sec.  II.  The  differences  in  the  eigen-structures  influence  the  formulation 
of  the  distribution  matrices  as  described  next.  Below,  we  first  describe  the  baseline  RD  scheme  for  the 
preconditioned  nonlinear  system,  and  then  upgrade  it  to  a  higher  order  by  the  correction  approach. 


B.  Baseline  RD  Scheme 


The  cell- residual  for  the  preconditioned  system  is  evaluated  in  two  steps.  First,  we  evaluate  the  cell-residuals 
of  the  nonlinear  equation  by  integrating  the  right  hand  side  of  Eq.  (91)  over  an  element  E  £  {E}  to  get  the 
unpreconditioned  cell- residuals 
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The  preconditioned  cell-residual,  which  is  distributed  to  the  element  vertices,  is  then  defined  as 


= 


^  U 

p 


=  PtyE, 


(97) 


where  P  is  evaluated  by  the  arithmetic  average  of  the  solution  U  within  the  element.  The  matrix  K j 
corresponding  to  the  preconditioned  system  is  defined  by 
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where  the  flux  Jacobians  are  evaluated  by  the  arithmetic  average  of  the  solution,  and  Kf  and  Kf  are, 
respectively,  the  contributions  from  the  advection  and  diffusion  components,  defined  as 
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and  v  is  the  arithmetic  average  of  the  diffusion  coefficient  over  the  element,  and 
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It  is  clear  that  the  eigenvalues  of  the  advective  and  diffusive  components  of  the  nonlinear  hyperbolic 
advection-diffusion  system  are  equivalent  to  the  linear  system: 


Aa  = 


\d  _ 

A1  -  “1 


(103) 


Similarly  to  the  linear  advection-diffusion  system,  we  construct  the  SUPG  distribution  matrix  of  the 
nonlinear  system  over  the  element  E  as: 


B  fUPG  =  Ii  +  D?  +  Df, 


(104) 


where 


and 


Df  = 


dfUPG  0  0 

0  0  0 
0  0  0 


dsupG  = 


D  d=1--^ 
1  2 v^3 


2  nrax(0,  ani)\ni\ 


,  Kf=-Aflni|, 


e?=ik r  2 


(105) 


(106) 


where  Kf+  is  constructed  based  on  the  projection  of  the  Af  onto  the  Af  running  wave.  We  note  here  that 
the  diffusive  flux  Jacobian  matrix  Af  has  the  following  eigenvectors: 


(107) 


R-n  = 

"  — 1/Ai 

—  1/ A2 

0 

T  —  1 

A1A2 
—  A1A2 

Al  fix 

A2  flx 

Al  hy 

—\2fly 

Tlx 

Tlx 

—  fly 

’  ~  A1-A2 

_  hy 

fly 

f^X 

0 

—  (Ai  —  A2)ny 

1 

(M 

1 
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which  is  different  from  the  linear  case.  Thus,  the  projection  matrix  II2  for  the  nonlinear  case  becomes 


n2 


Ai  —  A2 


Ai 

-A1A2  fix 
-A1A2  fly 


-A  2fll  -A  2  flx  fly 

-A  2nxhy  -A  2fiy 


and  the  K^+  can  be  expressed  as 


K?+  =  XllU,  n,  . 


(108) 


(109) 


The  residual  at  a  node  j  is  then  defined  in  the  form  of  Eq.  (10).  This  completes  the  construction  of  the 
baseline  RD  scheme  for  the  nonlinear  advection-diffusion  equation. 


C.  RD-CC2  and  RD-CC3  Schemes 

The  high-order  schemes,  RD-CC2  and  RD-CC3,  can  be  constructed  by  improving  the  cell-residuals  with  the 
high-order  curvature  correction  approach  discussed  in  Secs.  V  and  VI.  Special  attention  should  be  paid, 
however,  in  constructing  the  high-order  hyperbolic  RD  schemes  for  nonlinear  equations:  the  A p  and  A q 
terms  in  Eq.  (71)  require  the  solution  gradients,  while  p  and  q  variables  in  the  nonlinear  formulation  are  the 
diffusive  fluxes.  With  this  in  mind,  the  construction  is  straightforward. 

The  cell-residuals  of  the  RD-CC3  scheme  for  the  nonlinear  advection-diffusion  equation  are  given  by 


= 


-  \  E  [  {st  -  %)  +  (*?  -  w 


i—  1 


_  3 

V 


i= 1 
3 


=  *i+^r'E[5inv*  +  W''dnE 


where  Sf,  S f  and  <5?  are  given  by  Eqs.  (78),  (79),  and  (80),  respectively,  and 

1 

6 ? 


ST 


xpX 


rq/l 


g  (A/x.A^i  +  AfyiAyi) , 
g  {AyXiA T  AgyiAyi) , 

i  (A(p/v)iAxi  +  A{q/u)iAyi) , 

1 

{A{p/v)XiAxi  +  A(jp/v)ViAyi) , 
(A(q/u)XiAxi  +  A(q/v)yiAyi) . 


6 


(110) 

(111) 

(112) 

(113) 

(114) 

(115) 

(116) 
(117) 


For  the  flux  gradients  in  the  curvature  correction  terms  for  /  and  g ,  we  employ  the  chain  rule  to  avoid  the 
LSQ  gradient  reconstruction: 


f  df  df 

h  =  duU‘  =  du{rM’ 


f  9f  df 
1 »  =  duu’  =  On (qM’ 


(118) 


and  similarly  for  the  flux  g.  The  gradients  of  s,  p,  q ,  p/v ,  q/v,  and  the  source  term,  s,  are  computed  by 
the  quadratic  LSQ  fit.  The  RD-CC3  scheme  is  obtained  by  distributing  the  above  residuals  by  the  SUPG 
distribution  matrix.  The  RD-CC2  scheme  is  obtained  simply  by  removing  Sf ,  <5|,  Sf,  5PJV ,  and  SqJv  terms 
from  the  cell-residuals.  Consequently,  the  RD-CC2  scheme  does  not  require  any  LSQ  gradient  reconstruction 
and  therefore,  is  compact  as  in  the  linear  case.  We  confirmed  by  numerical  experiments  that  the  RD-CC2  and 
RD-CC3  schemes  preserve  exact  quadratic  and  cubic  solutions,  respectively,  by  the  method  of  manufactured 
solutions. 
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VIII.  Results 


In  this  section,  we  verify  and  examine  the  accuracy  of  the  hyperbolic  RD  schemes.  In  all  the  test  problems, 
we  specify  the  exact  solution  u  on  the  boundaries,  and  initialize  all  the  variables  in  the  interior  of  the  domain 
with  a  zero  value,  and  solve  the  hyperbolic  system  for  u ,  p ,  and  q.  The  linear  relaxation  is  performed  with  a 
Gauss-Seidel  (GS)  algorithm  to  reduce  the  linear  residual  by  five  orders  of  magnitude  with  the  maximum  of 
1000  relaxations.  Typically,  it  takes  300-600  relaxations  with  the  under-relaxation  parameter  in  the  range  of 
0.5  to  0.8.  The  implicit  solver  is  taken  to  be  converged  at  ten  orders  of  magnitude  residual  reduction  for  all 
equations;  this  is  obtained  typically  with  3-5  Newton  iterations  for  the  RD-CC2  scheme,  and  10-15  Newton 
iterations  for  the  RD-CC3  scheme. 

Simulations  were  performed  using  both  the  unified  and  the  non-unified  approaches  discussed  in  Sec.  III. 
The  convergence  speed  and  the  order  of  accuracy  results  were  found  to  be  nearly  identical  with  either  of  the 
two  approaches.  Therefore,  we  only  present  the  results  obtained  with  the  non-unified  approach. 

We  also  note  that  for  some  anisotropic,  stretched  grids,  we  could  not  obtain  solution  with  LDA  because 
the  GS  linear  solver  was  unstable.  A  more  efficient  and  stable  linear  solver  is  needed,  particularly  for  the 
LDA  scheme,  but  we  did  not  pursue  this.  Therefore,  we  only  focus  on  the  hyperbolic  RD  schemes  with 
SUPG  distribution  method. 

A.  Linear  Advection-Diffusion  Problem 

For  the  linear  advection-diffusion  equation  (1),  we  verify  the  order  of  accuracy  and  the  quality  of  the  predicted 
solution  and  solution  gradients  for  isotropic  and  anisotropic  irregular  triangular  grids,  and  for  a  problem 
with  curved  boundaries. 


1.  Isotropic  Grids 

Consider  a  steady  two-dimensional  linear  advection-diffusion  equation,  Eq.  (1),  with  s  =  0.  This  equation 
has  an  exact  exponential  solution28  of  the  form 


u(x,  y)  =  C  cos(Airri)  exp 

where  A  and  C  are  arbitrary  constants,  and  £  =  ax  +  by  and  p  =  bx  —  ay.  This  solution  is  smooth  and  has 
no  singularly  on  the  boundaries  and  therefore  is  appropriate  for  the  order  of  accuracy  verification. 

A  series  of  isotropic  irregular  grids  was  generated  and  the  accuracy  of  the  proposed  high-order  schemes  is 
verified  against  the  baseline  RD  and  RD-JCP2010  schemes.  The  first  three  coarse  grids  are  shown  in  Fig.  5. 
Error  convergence  results  are  shown  in  Fig.  6  along  with  a  more  detailed  information  in  Tables  5-7.  For 
the  variable  u,  all  second-order  schemes  give  truly  second-order  accuracy,  and  the  RD-CC3  scheme  gives 
fourth-order  accuracy.  For  the  gradient  variables,  p  and  q,  on  the  other  hand,  the  RD-CC2  and  RD-CC3 
schemes  yield,  respectively,  second-  and  third-order  accuracy,  but  the  baseline  RD  and  RD-JCP2010  do  not 
give  truly  second-order  accuracy. 


2v 


(119) 


Table  5.  LI  error  convergence  of  u  for  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular  grids 

(A  =  2,  C  =  -1,  a  =  2,  b  =  1,  v  =  0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

5.66E-05 

- 

6.98E-05 

- 

1.91E-05 

- 

2.31E-06 

- 

64x64 

1.08E-05 

2.39 

1.31E-05 

2.41 

4.11E-06 

2.22 

1.43E-07 

4.01 

128x28 

2.47E-06 

2.13 

2.82E-06 

2.22 

9.87E-07 

2.06 

1.02E-08 

3.81 

256x256 

5.73E-07 

2.11 

6.18E-07 

2.19 

2.38E-07 

2.02 

6.94E-10 

3.88 

Iterative  convergence  results  are  provided  in  Table  8  with  the  total  number  of  Newton  iterations  and  the 
number  of  GS  linear  relaxations  per  Newton  iteration.  As  expected,  the  implicit  solver,  which  is  Newton’s 
method  for  the  baseline,  the  RD-JCP2010  and  the  RD-CC2  schemes,  converges  with  3-5  Newton  iterations. 
For  the  RD-CC3  scheme,  the  implicit  solver,  which  is  not  truly  Newton’s  method,  converges  also  very  rapidly 
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(a)  8  x  8 


(b)  16  x  16 


(c)  32  x  32 


Figure  5.  Samples  of  the  irregular  grids  used  in  this  study. 


(a)  u 


(b)  p  =  ux 


(c)  q  =  uy 


Figure  6.  Order  of  accuracy  comparisons  of  the  baseline  RD,  RD-JCP2010,  RD-CC2,  and  RD-CC3  schemes. 


Table  6.  LI  error  convergence  of  p  (=  ux)  for  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular 
grids  (A  =  2,  C  =  —1,  a  =  2,  b  =  1,  is  =  0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

5.19E-03 

- 

7.87E-03 

- 

5.25E-04 

- 

7.43E-05 

- 

64x64 

1.70E-03 

1.61 

2.04E-03 

1.94 

1.20E-04 

2.13 

5.88E-06 

3.66 

128x128 

5.88E-04 

1.53 

5.54E-04 

1.88 

3.18E-05 

1.92 

5.85E-07 

3.33 

256x256 

1.99E-04 

1.56 

1.69E-04 

1.71 

9.14E-06 

1.80 

7.06E-08 

3.05 

Table  7.  LI  error  convergence  of  q  (=  uy)  for  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular 
grids  (A  =  2,  C  =  —1,  a  =  2,  b  =  1,  is  =  0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

7.04E-03 

- 

6.93E-03 

- 

6.16E-04 

- 

8.58E-05 

- 

64x64 

2.32E-03 

1.60 

2.20E-03 

1.65 

1.46E-04 

2.08 

6.61E-06 

3.70 

128x128 

7.67E-04 

1.60 

6.91E-04 

1.67 

3.82E-05 

1.93 

6.50E-07 

3.35 

256x256 

2.54E-04 

1.59 

2.14E-04 

1.69 

1.03E-05 

1.89 

7.62E-08 

3.09 
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Table  8.  Iterative  convergence  of  the  various  hyperbolic  RD  schemes  on  irregular  grids  (A  =  2,  C  =  —1,  a  =  2, 
6=l,i/  =  0.01).  Solutions  are  considered  converged  when  the  minimum  of  ten  orders  of  magnitude  reduction  is 
achieved  for  all  the  nodal-residuals.  The  Gauss-Seidel  relaxation  is  performed  for  each  Newton  iteration  until 
at  least  five  orders  of  magnitude  reduction  is  obtained  in  the  linear  residuals. 


Grids 

RD 

RD-JCP2010 

RD-CC2 

RD-CC3 

Newton  (Nt) 

GS/Nt 

Newton 

GS/Nt 

Newton 

GS/Nt 

Newton 

GS/Nt 

32x32 

3 

68 

3 

48 

3 

72 

10 

76 

64x64 

3 

103 

3 

69 

3 

76 

10 

83 

128x128 

3 

165 

3 

111 

3 

118 

10 

170 

256x256 

3 

273 

3 

190 

3 

249 

10 

361 

in  10  Newton  iterations.  Note  also  that  asymptotically  the  number  of  linear  relaxations  increases  linearly  (not 
quadratically  as  typically  in  conventional  methods)  with  the  grid  sizes.  This  is  a  feature  of  the  hyperbolic 
method,  arising  from  the  elimination  of  the  second  derivatives  in  the  target  equation.5 


(a)  ux\  RD  (baseline) 


(b)  ux:  RD-JCP2010 


(c)  ux  (exact) 


(d)  uy\  RD  (baseline) 


(e)  uy-,  RD-JCP2010 


(f)  uy  (exact) 


Figure  7.  Comparisons  between  solution  gradients  (ux  and  uy)  predicted  by  the  baseline  and  RD-JCP2010 
schemes  on  irregular  32x32  grid  (see  Fig.  5b)  with  the  exact  values  (A  =  2,  C  =  —0.009,  a  =  2,  6  =  1,  u  =  0.01). 


To  demonstrate  the  high  quality  of  the  predicted  solution  gradients  within  the  domain  by  the  improved 
RD-CC2  and  RD-CC3  schemes,  we  show  results  for  the  32x32  irregular  grid  shown  in  Fig.  5b.  Accurate 
prediction  of  the  solution  gradients  within  the  domain  is  of  great  interest  in  many  applications  (e.g.,  turbulent 
flows,  reacting  flows,  etc.).  Figure  7  shows  the  solution  gradients  predicted  by  the  baseline  RD  and  RD- 
JCP2010  schemes,  and  compares  them  with  the  exact  solutions.  As  can  be  seen  clearly,  the  predicted 
results  are  very  oscillatory  and  inaccurate.  On  the  other  hand,  as  shown  in  Fig.  8,  the  gradients  predicted 
by  the  proposed  RD-CC2  and  RD-CC3  schemes  are  very  accurate  and  smooth.  Note  that,  our  numerical 
experiments  suggest  that  in  some  cases  the  noise  produced  by  the  baseline  RD  and  RD-JCP2010  schemes 
will  not  completely  disappear  in  the  entire  domain  even  on  a  highly  refined  grid. 
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(a)  Ux'i  RD-CC2 


(b)  ux ;  RD-CC3 


(c)  ux  (exact) 


(d)  Uy ;  RD-CC2 


(e)  uy ;  RD-CC3 


(f)  uy  (exact) 


Figure  8.  Comparisons  between  solution  gradients  (ux  and  uy)  predicted  by  the  RD-CC2  and  RD-CC3  schemes 
on  irregular  32x32  grid  (see  Fig.  5b)  with  the  exact  values  ( A  =  2,  C  =  —0.009,  a  =  2,  6=1,  v  =  0.01). 
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The  reconstructed  gradients  computed  by  the  quadratic  LSQ  fits  from  the  solution  of  the  RD-CC3 
scheme,  compared  with  the  gradients  predicted  directly  with  the  RD-CC3  scheme  are  shown  in  Figs.  9  and 
10.  Clearly  the  reconstructed  gradients  are  less  accurate  and  deteriorate  even  more  within  the  domain.  This 
comparison  emphasizes  the  difficultly  in  obtaining  accurate  solution  gradients  by  reconstruction  techniques, 
a  common  practice  with  conventional  schemes. 


Figure  9.  Comparison  of  ux  reconstructed  with  the  quadratic  LSQ  fit  from  the  solution  of  the  RD-CC3  scheme 
and  the  one  predicted  directly  with  the  RD-CC3  scheme  on  irregular  32x32  grid  in  Fig.  5b.  (A  =  2,  C  =  —0.009, 
a  =  2,  b  =  1,  v  =  0.01). 


2.  Uniform  Accuracy  over  Re  =  a/v 

In  this  section  we  demonstrate  uniform  accuracy  predicted  by  the  RD-CC2  and  RD-CC3  schemes  for  both 
the  solution  and  the  solution  gradients  over  ranges  of  Re  with  different  irregular  grids.  These  results  are 
shown  in  Fig.  11,  which  verifies  that  the  RD-CC2  and  RD-CC3  schemes  produce,  respectively,  at  a  minimum, 
second-  and  third-order  solutions  for  all  the  variables.  Note  that  in  the  advection-limit  {v  — >  0)  the  hyperbolic 
advection-diffusion  system  reduces  to  a  simple  scalar  advection  equation,  for  which  the  RD-CC2  scheme  is 
third-order  accurate  for  all  the  variables.  In  the  diffusion-limit  ( v  -A  oo),  the  hyperbolic  advection-diffusion 
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(c)  Exact;  uy 


(d)  uy  at  y  =  0 


Figure  10.  Comparison  of  uy  reconstructed  with  the  quadratic  LSQ  fit  from  the  solution  of  the  RD- CC3 
scheme  and  the  one  predicted  directly  with  the  RD— CC3  scheme  on  irregular  32x32  grid  in  Fig.  5b.  (A  =  2, 

C  =  -0.009,  a  =  2,  b  =  1,  v  =  0.01). 
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RD-CC2,  Irregular  Grid 


RD-CC2,  Irregular  Grid 


(a)  RD-CC2;  u 

RD-CC3, 1  rregular  Grid 


RD-CC3,  Irregular  Grid 


RD-CC2, 1  rregular  Grid 


RD-CC3, 1  rregular  Grid 


(f)  RD-CC3;  uv 


Figure  11.  Uniform  accuracy  over  ranges  of  Re  with  the  proposed  RD-CC2  and  RD-CC3  schemes  on  irregular 
grids  (A  =  2,  C  =  —0.009,  a  =  2,  b  ==  1,  Re  = 


system  becomes  a  hyperbolic  diffusion  system,  for  which  the  RD-CC3  scheme  is  fourth-order  for  u  but 
remains  third-order  for  the  solution  gradients. 

Figure  12  shows  the  Newton  iteration  versus  Re  for  both  RD-CC2  and  RD-CC3  schemes.  It  is  clear  that 
the  number  of  Newton  iterations  is  independent  of  the  Re  and  grid  sizes.  We  remark  that  some  instabilities 
were  observed  with  the  linear  relaxation  for  high-Re  cases  ( Re  >  103) ,  but  found  that  both  RD-CC2  and 
RD-CC3  schemes  become  stable  if  the  curvature  correction  term  applied  to  the  advection  fluxes  (i.e.,  5f  ), 
is  evaluated  with  quadratic  LSQ  reconstructions;  that  is,  ux  and  uy  variables  present  in  the  <5“  term  (see 
Eq.  71)  are  reconstructed  from  the  solution  variable  u,  instead  of  being  substituted  by  the  gradient  variables 
p  and  q.  The  jump  in  the  number  of  Newton  iterations  shown  in  Fig.  12  for  the  RD-CC2  is  due  to  this 
switch,  which  is  implemented  to  avoid  instability  for  high-i?e  cases.  We  also  remark  that  if  the  switch  is 
activated  for  all  the  Re,  the  jump  will  disappear  and  a  similar  result  shown  for  the  RD-CC3  scheme  will  be 
recovered. 

A  linear  convergence,  O(N),  was  also  obtained  for  the  ranges  of  Re  and  is  shown  in  Fig.  13.  For  high-i?e 
cases,  a  better-than- linear  convergence  is  obtained  because  of  high  nresh-i?e  (i.e.,  ah/v)  values  used  with 
coarser  grids.  For  example,  for  the  v  =  10~3  case  (Re  =  1000\/5),  linear  convergence  is  obtained  from  the 
second  coarsest  grid  (mesh-i?e  <  35).  Note  that  the  RD-CC2  and  RD-CC3  schemes  are  efficient  enough  in 
converging  high-i?e  cases  on  grids  with  high  rnesh-i?,e,  and  therefore  we  did  not  attempt  to  fix  a  mesh-Re 
on  these  cases. 

3.  High  A  sped- Ratio  Grids 

In  practical  applications,  we  are  often  interested  in  the  gradients  at  the  boundary  and  with  highly  anisotropic 
grids,  e.g.,  viscous  drag  or  heat  flux  predictions  in  high-Reynolds-number  viscous  flows.  To  demonstrate  the 
benefit  of  the  hyperbolic  schemes  for  such  applications,  we  consider  two  highly  stretched  anisotropic  irregular 
grids  (see  Fig.  14).  The  y-derivative  of  u  at  y  =  0  is  plotted  and  compared  with  the  exact  solution  in  Fig.  14. 
The  predicted  normal  gradient  is  very  accurate  and  smooth  on  such  relatively  coarse  irregular  grids.  Similar 
results  have  been  obtained  by  other  schemes,  and  therefore  are  not  shown. 
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(a)  RD-CC2 


(b)  RD-CC3 


Figure  12.  Newton  iterations  used  by  the  proposed  RD-CC2  and  RD-CC3  schemes  on  irregular  grids  over 
ranges  of  Re  and  grids  (A  =  2,  C  =  —0.009,  a  =  2,  6  =  1,  Re  =  '/E/u). 


RD-CC3, 1  rregular  Grid 


'°9ioh 

(a)  RD-CC2 


Figure  13.  Average  Gauss-Seidel  relaxations  per  Newton  iteration  used  by  the  proposed  RD-CC2  and  RD-CC3 
schemes  on  irregular  grids  over  ranges  of  Re  and  grids  (A  =  2,  C  —  —0.009,  a  =  2,  6=1,  Re  =  v5/V). 
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Derivatives  of  u  at  y=0 
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(b)  32  x  32 


(c)  16  x  16 


(d)  32  x  32 


Figure  14.  Anisotropic,  stretched,  and  irregular  grids,  and  the  corresponding  predictions  of  the  normal 
gradients  along  the  y  =  0  boundary  with  the  baseline  hyperbolic  RD  scheme  (A  =  2,  C  =  —0.009,  a  =  2,  b  =  1, 
v  =  0.01). 
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This  exceptional  feature  of  the  hyperbolic  RD  schemes  is  further  illustrated  by  comparing  the  normal 
gradients  predicted  by  the  baseline  RD  scheme  and  a  conventional  FV  scheme  on  a  high  aspect-ratio  (AR) 
grid.  The  grid  and  the  results  are  shown  in  Fig.  15  where  the  normal  gradients  at  the  y  =  0  and  the  x  =  0 
boundaries  predicted  with  the  hyperbolic  RD  scheme  are  compared  with  those  obtained  by  the  linear  LSQ 
reconstruction  from  the  solution  of  a  conventional  FV  scheme  taken  from  Ref.  8,  where  the  scalar  advection- 
diffusion  equation  is  solved  by  the  second-order  edge-based  advection  scheme  and  the  second-order  Galerkin 
discretization  for  diffusion  (referred  to  as  Galerkin  in  the  figure).  The  predicted  normal  derivative  using  the 
baseline  RD  scheme  is  remarkably  smooth  and  accurate,  while  the  gradients  obtained  by  the  conventional 
scheme  shows  oscillatory  behavior  and  inaccurate  predictions.  For  this  problem  also,  similar  results  have 
been  obtained  by  other  schemes,  and  therefore  are  not  shown. 

These  results  indicate  that  even  the  baseline  RD  and  RD-JCP2010  schemes  can  produce  accurate  gra¬ 
dients  on  boundaries  although  they  yield  oscillations  in  the  interior  as  shown  in  the  previous  section.  Also, 
the  results  demonstrate  that  the  RD-CC2  and  RD-CC3  schemes  perform  successfully  even  on  irregular 
high-aspect-ratio  grids,  which  are  representatives  of  adapted  viscous  grids. 


X 


(a)  High  AR  grid  (b)  uy  (c)  ux 

Figure  15.  Comparisons  between  the  gradients  (ux  and  uy)  predicted  with  the  baseline  hyperbolic  RD  scheme 
and  the  published  result  of  the  conventional  Galerkin  method8  on  a  high  aspect-ratio  grid  with  AR  =  100; 
the  inset  figure  shows  a  portion  of  the  grid  with  1:1  scaling  on  both  axes  (A  =  2,  C  =  1,  a  =  1.23,  b  =  0.12, 
v  =  0.1235839795). 


B.  Domain  with  Curved  Geometrical  Boundaries 

In  this  section,  we  verify  the  accuracy  of  the  RD-CC2  and  RD-CC3  schemes  for  a  problem  with  a  curved 
boundary.  The  problem  is  a  potential  flow  over  a  unit  circle,  taken  from  Ref.  29.  The  governing  equation  is 
the  Laplace  equation  for  the  stream  function,  ip: 

dxxip  +  dyyip  =  0.  (120) 

We  solve  the  above  equation  by  solving  the  following  hyperbolic  system  of  equations: 

dTtp  =  dxp  +  dyq,  dTp  =  (dxip  -  p) ,  dTq  =  ( dyip  -  q) ,  (121) 

where  the  x-  and  y- velocity  components  are  related  to  the  gradient  of  the  main  variable  ip",  i.e.,  tpx  =  — 1> 
and  ipy  =  u. 

We  discretize  the  hyperbolic  potential  flow  system  of  equations  by  the  RD-CC2  and  RD-CC3  schemes, 
described  in  Secs.  V  and  VI,  and  verify  their  predicted  orders  of  accuracy  over  a  set  of  eight  irregular 
anisotropic  triangular  grids  of  441,  1681,  3721,  6561,  10201,  14641,  25921,  and  40401  nodes.  The  grids  are 
composed  of  linear  elements  only,  and  no  curved  elements  are  used.  The  domain  and  a  sample  grid  (1681 
nodes)  is  shown  in  Fig.  16a.  We  impose  the  exact  solution  on  the  outer  domain,  and  specify  ip  =  0  at 
the  solid  boundary  nodes  to  focus  on  the  order  of  accuracy  in  the  gradients.  Note  that  the  gradients  at 
boundary  nodes  are  predicted  by  the  numerical  schemes.  Figure  16  shows  the  contours  of  ipx  obtained  with 
the  RD-CC3  scheme  on  the  1681-node  grid.  The  excellent  quality  of  the  predicted  derivatives  is  observed  in 
comparison  with  the  exact  solution.  The  predicted  streamline  is  also  shown  in  the  same  figure. 
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(c)  Predicted  streamline  on  the  grid  of  1681  nodes.  (d)  Predicted  (—'i/’x)  on  the  grid  of  1681  nodes. 


Figure  16.  Potential  flow  (i/jXx  =  0)  over  a  unit  circle  on  irregular  and  stretched  grids  with  RD-CC3.  The 

contours  show  the  ^-gradient  of  the  solution  variable  i/>. 
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The  order  of  accuracy  of  theses  schemes  on  predicting  solution  gradients  on  curved  boundaries  is  also 
verified.  The  error  convergence  results  for  the  boundary  nodes  and  the  entire  domain  are  shown  in  Fig.  17. 
As  it  is  shown,  the  design  order  of  accuracy  is  clearly  observed  through  the  boundary  nodes.  It  is  also 

Potenti  al  f  I  ow,  I  rregul  ar  &  Stretched  G  ri  d  Potenti  al  f  I  ow,  I  rregul  ar  &  Stretched  G  ri  d 


Figure  17.  Potential  flow  (ipXx  +  ^yy  =  0)  over  a  unit  circle  on  several  irregular  and  stretched  grids  with  the 
RD-CC2  and  RD-CC3  schemes. 

evident  that  the  RD-CC3  scheme  is  more  accurate  than  the  RD-CC2  scheme,  even  on  the  coarsest  grid 
level,  indicating  that  the  third-order  scheme,  which  is  constructed  for  linear  elements,  is,  in  fact,  solving  the 
correct  problem  on  the  curved  boundary  that  is  defined  only  with  linear  elements. 

C.  Nonlinear  Advection-Diffusion  Problem 

Consider  the  following  nonlinear  advection-diffusion  equation: 

dtu  +  dxf  +  dyg  =  dx(ydxu )  +  dy{vdyu)  +  s(x,  y),  (122) 

where  /  =  g  =  u2 / 2,  v  =  u.  The  source  term  s(x,  y )  is  defined  by 

s(x,  y)  =  ue(uex  +  uey)  -  ( u% )2  -  ( uey )2  -  ue(uexx  +  ueyy),  (123) 


+  C0,  (124) 

where  C  =  —0.009,  A  =  2.0,  =  0.01,  and  Co  =  1,  so  that  ue  is  the  exact  solution  to  the  nonlinear 

advection-diffusion  equation  (122).  Note  that  Co  must  be  greater  than  C  in  order  for  the  diffusion  coefficient 
to  be  positive. 

We  solve  the  nonlinear  advection-diffusion  equation  in  a  square  domain  of  [0,1]  x  [0,1],  and  verify  the  order 
of  accuracy  of  the  solution  and  solution  gradients  predicted  by  the  baseline  RD,  the  RD-CC2,  and  the  RD- 
CC3  schemes.  The  order  of  accuracy  comparison  is  shown  in  Fig.  18  along  with  more  detailed  information  in 
Tables  9-11,  illustrating  that  the  RD-CC2  and  RD-CC3  schemes  achieve  an  improved  order  of  accuracy  for 
both  the  solution  and  the  solution  gradients  on  irregular  grids.  On  the  other  hand,  the  baseline  RD  scheme 
suffers  from  the  reduced  order  of  accuracy  in  the  gradients. 


=  C  cos(Anri)  exp 


1  —  \Jl  +  4A2Tr2is{ 


2vn 


£ 


IX.  Conclusions 

We  have  presented  detailed  formulation  and  implementation  procedures  for  hyperbolic  RD  schemes  for 
general  advection-diffusion  problems  on  arbitrary  triangular  grids.  Improving  upon  the  previous  work,5,6 
we  developed  new  second-  and  third-order  hyperbolic  RD  schemes  that  preserve,  respectively,  quadratic  and 
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Figure  18.  Order  of  accuracy  comparisons  between  the  baseline  and  the  proposed  hyperbolic  RD  schemes  on 
irregular  grids  for  the  nonlinear  advection-diffusion  problem  (f  =  g  =  u2  / 2,  v  =  u). 


Table  9.  L\  error  convergence  of  u  for  the  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular 
grids,  for  the  nonlinear  advection-diffusion  problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

7.93E-05 

- 

8.79E-05 

- 

3.98E-06 

- 

64x64 

1.93E-05 

2.04 

2.15E-05 

2.03 

2.43E-07 

4.04 

128x128 

4.90E-06 

1.98 

5.46E-06 

1.98 

1.59E-08 

3.93 

256x256 

1.22E-06 

2.01 

1.36E-06 

2.01 

1.12E-09 

3.83 

Table  10.  L\  error  convergence  of  ux  for  the  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular 
grids,  for  the  nonlinear  advection-diffusion  problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

1.00E-03 

- 

5.66E-04 

- 

5.84E-05 

- 

64x64 

3.63E-04 

1.46 

1.59E-04 

1.83 

5.94E-06 

3.30 

128x128 

1.44E-04 

1.33 

4.84E-05 

1.72 

7.75E-07 

2.94 

256x256 

5.31E-05 

1.44 

1.47E-05 

1.72 

1.12E-07 

2.79 

Table  11.  L\  error  convergence  of  uy  for  the  hyperbolic  RD  schemes  with  the  SUPG  distribution  on  irregular 
grids,  for  the  nonlinear  advection-diffusion  problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32x32 

1.46E-03 

- 

6.59E-04 

- 

6.17E-05 

- 

64x64 

5.62E-04 

1.38 

1.71E-04 

1.95 

6.28E-06 

3.30 

128x128 

2.23E-04 

1.33 

4.91E-05 

1.80 

7.94E-07 

2.98 

256x256 

8.23E-05 

1.44 

1.44E-05 

1.77 

1.14E-07 

2.80 
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cubic  solutions  on  arbitrary  triangular  grids.  We  also  developed  a  high-order  source  term  discretization, 
which  was  found  to  be  essential  in  constructing  the  third-order  scheme  with  linear  elements.  We  showed 
that  the  second-  and  third-order  schemes  can  be  constructed  by  a  separate  treatment  of  the  advective 
and  diffusive  terms,  enabling  development  of  hyperbolic  residual-distribution  schemes  for  the  compressible 
Navier-Stokes  equations.  The  superiority  of  the  proposed  schemes  was  further  demonstrated  with  accurate 
and  smooth  predictions  of  the  gradients  within  domains,  on  both  flat  and  curved  boundaries.  The  new  second- 
order  scheme,  which  is  still  defined  within  a  compact  stencil,  allows  an  efficient  construction  of  Newton’s 
method  with  the  exact  residual  Jacobian  computed  by  the  automatic  differentiation.  We  demonstrated  that 
the  resulting  implicit  solver  is  capable  of  reducing  the  residuals  of  all  equations  by  at  least  ten  orders  of 
magnitude  with  typically  less  than  10  Newton  iterations  for  the  second-order  scheme.  For  the  third-order 
scheme,  we  showed  that  the  linear  solver  converges  also  very  rapidly  with  10-15  Newton  iterations  with 
the  second-order  Jacobian.  Because  of  the  hyperbolic  reformulation  of  the  advect ion- diffusion  equation,  the 
number  of  linear  relaxations  per  Newton  iteration  was  also  shown  to  increase  linearly,  not  quadratically  as 
for  typical  diffusion  problems,  with  grid  refinement. 
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A.  Fourier  analysis 


In  this  appendix,  we  present  the  Fourier  analysis  to  investigate  the  stability  of  relaxation  schemes  applied 
to  the  linear  system  (14).  The  analysis  is  performed  for  the  baseline  RD  scheme  only,  but  the  results  are 
suggestive  and,  in  fact,  effective  to  stabilize  other  schemes. 

Consider  a  two-dimensional  stencil  shown  in  Fig.  19,  where  we  seek  the  solution  at  node  j,  which  is 


Figure  19.  Schematic  of  the  stencil  used  in  Fourier  analysis. 


surrounded  by  six  elements  and  and  six  nodes.  The  grid  aspect  radio  (AR)  is  defined  as  hx/hy.  We  now 
define  the  vectors  U  and  Q  to  all  the  nodes;  that  is 


u,= 


Pj 

<b 


Ufc  = 


uk 

Pk 

qk 


for  k  =  1..6, 


(125) 


Q?  _ 


o 

Pj  /Tr 
Qj/Tr 


Q  k  — 


0 

Pk/Tr 

qk/Tr 


for  k  =  1..6. 


The  cell-residual  then  becomes: 


3>t  =  —  ^2  KjUj  +  QT ST . 


i£T 


For  example,  for  the  element  containing  nodes  (j.  1,  2),  we  have: 

$Te(j’i’2)  =  -KjVj  -  K1XJ1  -  K2U2  +  ^(Qj  +  Qr  +  Q2)^yA 


(126) 


(127) 


(128) 


32  of  40 


American  Institute  of  Aeronautics  and  Astronautics 


We  remark  that  K  is  an  element  dependent  matrix  and  therefore  special  care  must  be  taken  in  evaluating 
the  elemental  residual  $T. 

Using  Eq.  (11),  we  evaluate  the  residual  of  the  node  j,  noting  that  the  median  dual  volume  Sj  for  our 
stencil  is  simply  hxhy.  We  now  start  the  Fourier  analysis  by  writing  the  Resj  as 


where  C j  and  C  k  are 


6 


;  =  CjVj  +  CfcU k, 

(129) 

k= 1 

<9Resj 

p  ^  j 

-  au,  ■ 

(130) 

9Res, 

( ~ i  J 

k  dVk  ' 

(131) 

We  then  transform  the  residual  ReSj  into  the  corresponding  Fourier  modes  by  replacing  the  vector  Ufc  with 
a  Fourier  mode  of  the  phase  angle  /?  =  (0X,  (3y)  with  j3  €  [0, 7r] : 


Ur  =  U,e 


i(/3x 


L+A 


Vk-Vj 
h‘ u 


(132) 


where  i  =  y/—l,  and  ( Xj,yj )  and  (xk,yk)  are  the  nodal  coordinates.  Equation  129  can  now  be  written  as: 

Res^  =  Uj  (Cj  +  e~i^C1  +  e~i^+^C2  +  e“i/3“C3  +  ei/3*C4  +  ei(~^+^C5  +  ei/J»C6)  (133) 
The  solution  of  Uj  using  the  Jacobi  relaxation  is  therefore  become: 

U"  • 1  =  MjBU"  (134) 

where  the  mass  matrix  M  is  defined  as 

Mjb  =  C  t1  (e-^Cr  +  e-i(-^+^)C2  +  e~^Cz  +  e^C4  +  e^+^Cs  +  ei/S»C6)  (135) 


Note  that  this  is  equivalent  to  the  Jabobi  relaxation  applied  to  the  linear  system  in  Eq.  (14)  without  the 
right  hand  side,  which  is  irrelevant  to  the  stability  analysis. 

The  three  complex  eigenvalues  of  the  matrix,  Mj^,  are  computed  for  several  RD  schemes  with  grid  AR  of 
1  and  106.  Here  we  used  50  equally  distributed  points  in  f3x  and  (3y  directions,  for  a  total  of  2500  points  (i.e., 
50x50).  We  first  present  the  results  of  the  first-order  Lax-Friedrichs  (or  Rusanov)  scheme30  (Fig.  20)  for  the 
hyperbolic  diffusion  problem.  The  x -  and  y- axes  are  the  real  and  imaginary  part  of  the  matrix  eigenvalues. 
The  colors  in  the  plots  correspond  to  the  three  eigenvalues.  This  Rusanov  scheme  is  stable  with  GS  for  any 
regular  grid.  Note  that  the  imaginary  component  of  the  eigenvalues  are  very  small. 

Figure  21  shows  the  results  for  the  hyperbolic  diffusion  system  using  the  LDA  scheme.  Clearly,  the 
results  exhibit  modes  with  magnitudes  larger  than  one,  specially  on  the  real  axis.  These  modes  indicate 
that  the  errors  could  amplify  and  may  not  be  damped.  Therefore,  the  LDA  scheme  is  unstable  with  the 
Jacobi  relaxation  for  our  system  equation.  Similar  results  are  obtained  with  the  Gauss-Seidel  relaxation. 
The  results  are  consistent  with  our  developed  two-dimensional  RD  code.  We  remark  that  while  the  figures 
are  shown  for  AR  =  106  are  identical  to  results  obtained  on  uniform  grids  with  arbitrary  AR.  Similar  results 
are  also  obtained  with  a  pure  hyperbolic  problem,  the  wave  equation,  utt  =  c2(uxx  +  uyy )  (see  Fig.  22).  The 
results  are  nearly  identical  to  the  hyperbolic  diffusion  problem  shown  in  Fig.  21  particularly  with  AR  =  1. 
The  results  indicate  that  the  instability  of  the  LDA  scheme  is  not  caused  by  the  first-order  hyperbolic  system 
method.  We  have  also  examined  the  stability  of  the  LDA  scheme  for  a  scalar  linear  advection  equation.  Note 
that  there  is  only  one  eigenvalue  for  the  pure  advection  problem.  These  results  are  presented  in  Fig.  23, 
which  indicates  non-damping  modes  specially  for  grids  with  AR  =  1.  The  results  confirm  that  the  presented 
instability  is  not  due  to  the  formulation  of  the  diffusion  problems  using  hyperbolic  first-order  system. 

The  SUPG,  LW,  and  stabilized-LDA  schemes  also  exhibit  similar  instability  with  the  GS  relaxation  for 
both  diffusion  and  scalar  linear  advection  problems  (see  Figs.  24-28).  Note  that  the  stabilized-LDA  scheme 
consists  of  the  LDA  scheme  and  a  stabilization  term  defined  as  K. 
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0.04 


-0.04  J 

(a)  Rusanov:  AR  =  1  (Diffusion)  (b)  Rusanov:  AR  =  106  (Diffusion) 

Figure  20.  Fourier  analysis  of  the  Rusanov  scheme  for  the  hyperbolic  diffusion  problem:  v  =  1,  hx  =  0.1. 


**  * 


*  * 


*****  * 

******%  * 

******  ** 
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(a)  LDA:  AR  =  1  (Diffusion) 


(b)  LDA:  AR  =  106  (Diffusion) 


Figure  21.  Fourier  analysis  of  the  LDA  scheme  for  the  hyperbolic  diffusion  system:  v  =  1,  hx  =  0.1. 
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(a)  LDA:  Af?,  =  1  (Wave) 


(b)  LDA:  A.R  =  106  (Wave) 


Figure  22.  Fourier  analysis  of  the  LDA  scheme  for  the  hyperbolic  wave  problem:  utt  =  c2(uxx  +uyy),  c  =  1, 
hx  =  0.1. 
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(a)  LDA:  AR  =  1  (Advection) 


(b)  LDA:  AR  =  10®  (Advection) 


Figure  23.  Fourier  analysis  of  the  LDA  scheme  for  linear  advection  equation:  a  =  2,  b  =  1,  hx  =  0.1. 
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(a)  SUPG:  AR  =  1  (Diffusion) 


(b)  SUPG:  AR  =  106  (Diffusion) 


Figure  24.  Fourier  analysis  of  the  SUPG  scheme  for  the  hyperbolic  diffusion  problem:  t/  =  1,  hx  =  0.1. 


(a)  SUPG:  Ai?  =  1  (Advection) 


(b)  SUPG:  AR  =  106  (Advection) 


Figure  25.  Fourier  analysis  of  the  SUPG  scheme  for  linear  advection  equation:  a  =  2,  b  =  1,  hx  =  0.1. 
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Figure  26.  Fourier  analysis  of  the  Lax-Wendrof  (LW)  scheme  for  the  hyperbolic  diffusion  problem:  u  =  1 
hx  =  0.1. 


-1-  -1- 
(a)  LW:  AR  =  1  (Advection)  (b)  LW:  AR  =  106  (Advection) 

Figure  27.  Fourier  analysis  of  the  Lax-Wendrof  (LW)  scheme  for  linear  advection  equation:  a  = 
hx  =  0.1. 
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(a)  Stabilized-LDA:  AR  =  1  (Diffusion)  (b)  Stabilized-LDA:  AR  =  10®  (Diffusion) 

Figure  28.  Fourier  analysis  of  the  stabilized-LDA  scheme  for  the  hyperbolic  advection-diffusion  problem:  u  =  1, 
hx  =  0.1. 


We  further  examined  the  instability  of  these  schemes  with  weighted  Jacobi  and  under-relaxed  GS  relax¬ 
ations.  The  matrix  for  the  weighted  Jacobi  relaxation,  TsAwjb  is 

Mivjb  =  (1-w)I  +  uMjb,  (136) 

where  I  is  an  identity  matrix  and  M^b  is  defined  in  Eq.  (135).  We  found  that  generally  these  schemes 
are  stable  with  weighted  Jacobi  or  under-relaxed  GS  relaxations  with  w  =  0.7  —  0.8.  For  instance,  Fig.  29 
presents  the  stability  plots  for  the  linear  advection  problem  using  SUPG  formulation.  The  results  indicate 
that  the  linear  solver  is  stable  on  arbitrary  regular  grids  with  oj  =  0.7  as  the  under-relaxation  value.  The  same 
technique  applied  to  the  hyperbolic  diffusion  problem  as  shown  in  Fig.  30.  Similar  results  were  obtained  with 
other  presented  schemes  and,  therefore,  not  repeated  here.  Note  that  in  the  Fourier  analysis,  effects  of  the 
boundary  points  and  grid  irregularity  could  not  be  included.  However,  we  found  that  the  under-relaxation 
parameter  of  uj  =  0.6  —  0.8  is  typically  sufficient  to  make  the  linear  relaxation  stable  on  arbitrary  anisotropic 
irregular  grids  for  general  advection-diffusion  problems. 
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(b)  SUPG:  AR  =  106,  cj  =  0.7  (Advection) 


Figure  29.  Fourier  analysis  of  the  SUPG  scheme  for  the  linear  advection  problem  (a  =  2,  b  =  1,  hx  =  0.1)  using 
weighted  Jacobi  iteration  (cj  =  0.7). 


Figure  30.  Fourier  analysis  of  the  SUPG  scheme  for  the  linear  advection  problem  ( v  =  1,  hx 
weighted  Jacobi  iteration  (cj  =  0.8). 
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